C foF2-hmF2-IRI.for
C-------------------------------------------------------------------- 
C shows how to use the subroutine iri_fpeak; note the COMMON
C blocks required
C
C irir_fpeak reads the following files:
C		IG_RZ.DAT indices file on unit=12
C		APF107.DAT indices files on unit=13 (ACCESS='DIRECT')
C       CCIR%%.ASC and URSI%%.ASC files on unit=10 (%%=month+10)
C		I/DGRF%%s.ASC on unit=14 (%%=year in 5-year increments)
C-------------------------------------------------------------------- 

        REAL		LATI,LONGI
      	COMMON 	/CONST/UMR  /const1/humr,dumr   /ARGEXP/ARGMAX
     &			/iounit/konsol,mess

		mess=.true.
		konsol=6
        ARGMAX=88.0
        pi=ATAN(1.0)*4.
        UMR=pi/180.
        humr=pi/12.
        dumr=pi/182.5

		print*,'year(yyyy),month(mm),day(dd),decimal UT hour'
		read(5,*) iyear,month,iday,hourut
		print*,'geographic latitude, longitude (0-360)'
		read(5,*) LATI,LONGI
        print*,'URSI(0) or CCIR(1)','storm off(0) on(1)'
        read(5,*) imap,istorm

		call iri_fpeak(iyear,month,iday,hourut,lati,longi,
     &		imap,istorm,foF2,foE,hmF2,xm3000)	

		print*,iyear,month,iday,hourut,lati,longi,
     &		foF2,xm3000,foE,hmF2
     	stop
     	end           


		subroutine iri_fpeak(iyear,month,iday,hourut,lati,longi,
     &		imap,istorm,foF2o,foEo,hmF2o,xm3000F2)	
C-------------------------------------------------------------------- 
C 
C INPUT:	iyear	year in yyyy format
C			month	month in mm format
C			day		day of month in dd format
C			hourut	Universal Time in decimal hours
C			lati	geographic latitude in degrees 
C 			longi	geographic longitude in degrees(0-360)
C			imap	=0 URSI, =1 CCIR maps for foF2
C			istorm	=0/1 storm model off/on
C
C OUTPUT:	foF2o	F2 peak plasma frequency in m-3
C			foEo	E peak plasma frequency in m-3
C			hmF2	F2 peak height in km
C			xm3000F2 M(3000)F2 factor 
C
C Reads IG_RZ.DAT indices file on unit=12
C		APF107.DAT indices files on unit=13 (ACCESS='DIRECT')
C       CCIR%%.ASC and URSI%%.ASC files on unit=10 (%%=month+10)
C		I/DGRF%%s.ASC on unit=14 (%%=year in 5-year increments)
C
C change to IRI-2014: COV for foE now F107_81
C-------------------------------------------------------------------- 
      
        REAL		LATI,LONGI,MODIP,MAGBR
      	CHARACTER	FILNAM*12
		DIMENSION	F2(13,76,2),FM3(9,49,2),F2N(13,76,2),FM3N(9,49,2),
     &				FF0(988),FF0N(988),XM0(441),XM0N(441),ARIG(3),
     & 				RZAR(3),indap(13)
     	LOGICAL		mess
      	COMMON 		/iounit/konsol,mess
               
		IUCCIR=10

        call MODA(0,iyear,MONTH,IDAY,IDAYNR,nrdaym)
        call tcon(iyear,month,iday,idaynr,rzar,arig,ttt,nmonth)
        if(nmonth.lt.0) goto 3330		! jump to end of program
        call APF_ONLY(iyear,month,iday,F107_daily,F107PD,F107_81,
     &        F107_365,IAP_daily)
        if(F107_daily.lt.0) goto 3330 
        cov=f107_81
c            f107y=f107PD
c            f10781=f107_81
c            f107365=f107_365
c			endif

        idayy=365
c  leap year rule: years evenly divisible by 4 are leap years, except
c  years also evenly divisible by 100 are not leap years, except years 
c  also evenly divisible by 400 are leap years. The year 2000 is a 100 
c  and 400 year exception and therefore it is a normal leap year. 
c  The next 100 year exception will be in the year 2100!
        if(iyear/4*4.eq.iyear) idayy=366    ! leap year
        ryear = iyear + (idaynr-1.0)/idayy
        call igrf_dip(lati,longi,ryear,300.0,dec,dip,magbr,modip)

c
c the program expects the coefficients files in ASCII format; if you
C want to use the binary version of the coefficients, please use the
C the statements that are commented-out below and comment-out the
C ASCII-related statements.
c
        	WRITE(FILNAM,104) MONTH+10
104         FORMAT('ccir',I2,'.asc')
        	OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &          FORM='FORMATTED')
        	READ(IUCCIR,4689) F2,FM3
4689    	FORMAT(1X,4E15.8)
        	CLOSE(IUCCIR)
c
c READ CCIR COEFFICIENT SET FOR NMONTH, i.e. previous 
c month if day is less than 15 and following month otherwise 
c
        	WRITE(FILNAM,104) NMONTH+10
        	OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &          FORM='FORMATTED')
        	READ(IUCCIR,4689) F2N,FM3N
        	CLOSE(IUCCIR)
c
c if URSI is selected than read URSI foF2 coefficients
c
        if(imap.lt.1) then
          	WRITE(FILNAM,1144) MONTH+10
1144        FORMAT('ursi',I2,'.asc')
          	OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &         FORM='FORMATTED')
          	READ(IUCCIR,4689) F2
          	CLOSE(IUCCIR)
          	WRITE(FILNAM,1144) NMONTH+10
          	OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &         FORM='FORMATTED')
          	READ(IUCCIR,4689) F2N
          	CLOSE(IUCCIR)
        endif
C
C LINEAR INTERPOLATION IN SOLAR ACTIVITY. IG12 used for foF2
C
        RR2=ARIG(1)/100.
        RR2N=ARIG(2)/100.
        RR1=1.-RR2
        RR1N=1.-RR2N
        DO 20 I=1,76
        DO 20 J=1,13
              K=J+13*(I-1)
              FF0N(K)=F2N(J,I,1)*RR1N+F2N(J,I,2)*RR2N
20            FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2

        RR2=RZAR(1)/100.
        RR2N=RZAR(2)/100.
        RR1=1.-RR2
        RR1N=1.-RR2N
        DO 30 I=1,49
        DO 30 J=1,9
            K=J+9*(I-1)
            XM0N(K)=FM3N(J,I,1)*RR1N+FM3N(J,I,2)*RR2N
30          XM0(K)=FM3(J,I,1)*RR1+FM3(J,I,2)*RR2

4292    zfof2  =  FOUT(MODIP,LATI,LONGI,HOURUT,FF0)
        fof2n  =  FOUT(MODIP,LATI,LONGI,HOURUT,FF0N)
        zm3000 = XMOUT(MODIP,LATI,LONGI,HOURUT,XM0)
        xm300n = XMOUT(MODIP,LATI,LONGI,HOURUT,XM0N)
        midm=15
        if(month.eq.2) midm=14
        if (iday.lt.midm) then
            fof2 = fof2n + ttt * (zfof2-fof2n)
            xm3000= xm300n+ ttt * (zm3000-xm300n)
        else
            fof2 = zfof2 + ttt * (fof2n-zfof2)
            xm3000= zm3000+ ttt * (xm300n-zm3000)
        endif			
        ABSLAT=ABS(LATI)
        hour=hourut+longi/15.	 		
        if(hour.gt.24.) hour=hour-24.
        CALL SOCO(idaynr,HOUR,LATI,LONGI,110.,SD,XHI,SAX,SUX)
        CALL SOCO(idaynr,12.0,LATI,LONGI,110.,SD,XHIN,SAX,SUX)
		FOE=FOEEDI(COV,XHI,XHIN,ABSLAT)
		rssn=rzar(3)
        HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000)
c
c stormtime updating for foF2 (foF2s, NmF2s) 
c
        if(istorm.gt.0) then
           	call apf(iyear,month,iday,hourut,indap)
            if(indap(1).gt.-1) then
            	icoord=1
            	kut=int(hourut)
        		call STORM(indap,lati,longi,icoord,
     &              cglat,kut,idaynr,stormcorr)
        		fof2=fof2*stormcorr     
        	endif
        endif

c		print*,iyear,month,iday,idaynr,hourut,lati,longi,
c     &		modip,foF2,xm3000,foE,hmF2           
		foF2o=foF2
		xm3000F2=xm3000
		foEo=FOE
		hmF2o=HMF2
		goto 3330
		
8448    if(mess) WRITE(konsol,8449) FILNAM
8449    FORMAT(1X////,
     &    ' The file ',A30,'is not in your directory.')

3330	CONTINUE
		RETURN
		END



C                     
C*************************************************************                  
C************* SUBROUTINES AND FUNCTIONS *********************                  
C*************************************************************                  
C
C
      real function FOUT(XMODIP,XLATI,XLONGI,UT,FF0)
c--------------------------------------------------------------
C CALCULATES CRITICAL FREQUENCY FOF2/MHZ USING SUBROUTINE GAMMA1.      
C XMODIP = MODIFIED DIP LATITUDE, XLATI = GEOG. LATITUDE, XLONGI=
C LONGITUDE (ALL IN DEG.), MONTH = MONTH, UT =  UNIVERSAL TIME 
C (DEC. HOURS), FF0 = ARRAY WITH RZ12-ADJUSTED CCIR/URSI COEFF.
C D.BILITZA,JULY 85.
c--------------------------------------------------------------
      DIMENSION FF0(988)
      INTEGER QF(9)
      DATA QF/11,11,8,4,1,0,0,0,0/
      FOUT=GAMMA1(XMODIP,XLATI,XLONGI,UT,6,QF,9,76,13,988,FF0)
      RETURN
      END
C
C
      real function XMOUT(XMODIP,XLATI,XLONGI,UT,XM0)
c--------------------------------------------------------------
C CALCULATES PROPAGATION FACTOR M3000 USING THE SUBROUTINE GAMMA1.
C XMODIP = MODIFIED DIP LATITUDE, XLATI = GEOG. LATITUDE, XLONGI=
C LONGITUDE (ALL IN DEG.), MONTH = MONTH, UT =  UNIVERSAL TIME 
C (DEC. HOURS), XM0 = ARRAY WITH RZ12-ADJUSTED CCIR/URSI COEFF.
C D.BILITZA,JULY 85.
c--------------------------------------------------------------
      DIMENSION XM0(441)
      INTEGER QM(7)
      DATA QM/6,7,5,2,1,0,0/
      XMOUT=GAMMA1(XMODIP,XLATI,XLONGI,UT,4,QM,7,49,9,441,XM0)
      RETURN
      END
C
C
      REAL FUNCTION HMF2ED(XMAGBR,R,X,XM3)         
c--------------------------------------------------------------
C CALCULATES THE PEAK HEIGHT HMF2/KM FOR THE MAGNETIC                           
C LATITUDE XMAGBR/DEGREE AND THE SMOOTHED ZUERICH SUNSPOT                         
C NUMBER R USING CCIR-M3000 XM3 AND THE RATIO X=FOF2/FOE.
C FOLLOWING CCIR RECOMMENDATION X IS LIMITED TO VALUE
C GREATER OR EQUAL TO 1.7 .                       
C [REF. D.BILITZA ET AL., TELECOMM.J., 46, 549-553, 1979]                       
C D.BILITZA,1980.     
c--------------------------------------------------------------
      F1=0.00232*R+0.222                         
      F2=1.2-0.0116*EXP(0.0239*R)            
      F3=0.096*(R-25.0)/150.0                      
      F4=1.0-R/150.0*EXP(-XMAGBR*XMAGBR/1600.0)
      if(x.lt.1.7) x=1.7
      DELM=F1*F4/(X-F2)+F3                
      HMF2ED=1490.0/(XM3+DELM)-176.0 
      RETURN          
      END             
C
C
C
        REAL FUNCTION GAMMA1(SMODIP,SLAT,SLONG,HOUR,
     &                          IHARM,NQ,K1,M,MM,M3,SFE)      
C---------------------------------------------------------------
C CALCULATES GAMMA1=FOF2 OR M3000 USING CCIR NUMERICAL MAP                      
C COEFFICIENTS SFE(M3) FOR MODIFIED DIP LATITUDE (SMODIP/DEG)
C GEOGRAPHIC LATITUDE (SLAT/DEG) AND LONGITUDE (SLONG/DEG)  
C AND UNIVERSIAL TIME (HOUR/DECIMAL HOURS). IHARM IS THE MAXIMUM
C NUMBER OF HARMONICS USED FOR DESCRIBING DIURNAL VARIATION.
C NQ(K1) IS AN INTEGER ARRAY GIVING THE HIGHEST DEGREES IN 
C LATITUDE FOR EACH LONGITUDE HARMONIC WHERE K1 GIVES THE NUMBER 
C OF LONGITUDE HARMONICS. M IS THE NUMBER OF COEFFICIENTS FOR 
C DESCRIBING VARIATIONS WITH SMODIP, SLAT, AND SLONG. MM IS THE
C NUMBER OF COEFFICIENTS FOR THE FOURIER TIME SERIES DESCRIBING
C VARIATIONS WITH UT.
C M=1+NQ(1)+2*[NQ(2)+1]+2*[NQ(3)+1]+... , MM=2*IHARM+1, M3=M*MM  
C SHEIKH,4.3.77.      
C---------------------------------------------------------------
      REAL*8 C(12),S(12),COEF(100),SUM             
      DIMENSION NQ(K1),XSINX(13),SFE(M3)           
      COMMON/CONST/UMR
      HOU=(15.0*HOUR-180.0)*UMR                    
      S(1)=SIN(HOU)   
      C(1)=COS(HOU)   

      DO 250 I=2,IHARM                             
        C(I)=C(1)*C(I-1)-S(1)*S(I-1)                 
        S(I)=C(1)*S(I-1)+S(1)*C(I-1)                 
250     CONTINUE        

      DO 300 I=1,M    
        MI=(I-1)*MM     
        COEF(I)=SFE(MI+1)                            
        DO 300 J=1,IHARM                             
          COEF(I)=COEF(I)+SFE(MI+2*J)*S(J)+SFE(MI+2*J+1)*C(J)                       
300       CONTINUE        

      SUM=COEF(1)     
      SS=SIN(SMODIP*UMR)                           
      S3=SS           
      XSINX(1)=1.0    
      INDEX=NQ(1)     

      DO 350 J=1,INDEX                             
        SUM=SUM+COEF(1+J)*SS                         
        XSINX(J+1)=SS   
        SS=SS*S3        
350     CONTINUE        

      XSINX(NQ(1)+2)=SS                            
      NP=NQ(1)+1      
      SS=COS(SLAT*UMR)                             
      S3=SS           

      DO 400 J=2,K1   
        S0=SLONG*(J-1.)*UMR                          
        S1=COS(S0)      
        S2=SIN(S0)      
        INDEX=NQ(J)+1   
        DO 450 L=1,INDEX                             
          NP=NP+1         
          SUM=SUM+COEF(NP)*XSINX(L)*SS*S1              
          NP=NP+1         
          SUM=SUM+COEF(NP)*XSINX(L)*SS*S2              
450       CONTINUE        
        SS=SS*S3        
400     CONTINUE
        
      GAMMA1=SUM      

      RETURN          
      END 
C
C
        REAL FUNCTION FOEEDI(COV,XHI,XHIM,XLATI)
C-------------------------------------------------------
C CALCULATES FOE/MHZ BY THE EDINBURGH-METHOD.      
C INPUT: 
C 	COV		MONTHLY MEAN 10.7CM SOLAR RADIO FLUX measured at 
C           ground level  
C   XHI		SOLAR ZENITH ANGLE IN DEGREE 
C   XHIM	SOLAR ZENITH ANGLE AT NOON IN DEGREE
C   XLATI 	ABSOLUTE VALUE OF GEOGRAPHIC LATITUDE IN DEGREE, 
C REFERENCE: 
C       KOURIS-MUGGELETON, CCIR DOC. 6/3/07, 1973
C       TROST, J. GEOPHYS. RES. 84, 2736, 1979 (was used
C               to improve the nighttime varition)
C       RAWER AND BILITZA, Adv. Space Res. 10(8), 5-14, 1990
C D.BILITZA--------------------------------- AUGUST 1986.    
        COMMON/CONST/UMR
C variation with solar activity (factor A) ...............
        A=1.0+0.0094*(COV-66.0)                      
C variation with noon solar zenith angle (B) and with latitude (C)
        SL=COS(XLATI*UMR)
        IF(XLATI.LT.32.0) THEN
                SM=-1.93+1.92*SL                             
                C=23.0+116.0*SL                              
        ELSE
                SM=0.11-0.49*SL                              
                C=92.0+35.0*SL  
        ENDIF
        if(XHIM.ge.90.) XHIM=89.999
        B = COS(XHIM*UMR) ** SM
C variation with solar zenith angle (D) ..........................        
        IF(XLATI.GT.12.0) THEN
                SP=1.2
        ELSE
                SP=1.31         
        ENDIF
C adjusted solar zenith angle during nighttime (XHIC) .............
        XHIC=XHI-3.*ALOG(1.+EXP((XHI-89.98)/3.))   
        D=COS(XHIC*UMR)**SP       
C determine foE**4 ................................................
        R4FOE=A*B*C*D     
C minimum allowable foE (foe_min=sqrt[SMIN])...............................
        SMIN=0.121+0.0015*(COV-60.)
        SMIN=SMIN*SMIN
        IF(R4FOE.LT.SMIN) R4FOE=SMIN                     
        FOEEDI=R4FOE**0.25                           
        RETURN          
        END   
C
C
        subroutine soco (ld,t,flat,Elon,height,
     &          DECLIN, ZENITH, SUNRSE, SUNSET)
c--------------------------------------------------------------------
c       s/r to calculate the solar declination, zenith angle, and
c       sunrise & sunset times  - based on Newbern Smith's algorithm
c       [leo mcnamara, 1-sep-86, last modified 16-jun-87]
c       {dieter bilitza, 30-oct-89, modified for IRI application}
c
c in:   ld      local day of year
c       t       local hour (decimal)
c       flat    northern latitude in degrees
c       elon    east longitude in degrees
c		height	height in km
c
c out:  declin      declination of the sun in degrees
c       zenith      zenith angle of the sun in degrees
c       sunrse      local time of sunrise in hours 
c       sunset      local time of sunset in hours 
c-------------------------------------------------------------------
c
        common/const/   dtr     /const1/humr,dumr
c amplitudes of Fourier coefficients  --  1955 epoch.................
        data    p1,p2,p3,p4,p6 /
     &  0.017203534,0.034407068,0.051610602,0.068814136,0.103221204 /
c
c s/r is formulated in terms of WEST longitude.......................
        wlon = 360. - Elon
c
c time of equinox for 1980...........................................
        td = ld + (t + Wlon/15.) / 24.
        te = td + 0.9369
c
c declination of the sun..............................................
        dcl = 23.256 * sin(p1*(te-82.242)) + 0.381 * sin(p2*(te-44.855))
     &      + 0.167 * sin(p3*(te-23.355)) - 0.013 * sin(p4*(te+11.97))
     &      + 0.011 * sin(p6*(te-10.41)) + 0.339137
        DECLIN = dcl
        dc = dcl * dtr
c
c the equation of time................................................
        tf = te - 0.5
        eqt = -7.38*sin(p1*(tf-4.)) - 9.87*sin(p2*(tf+9.))
     &      + 0.27*sin(p3*(tf-53.)) - 0.2*cos(p4*(tf-17.))
        et = eqt * dtr / 4.
c
        fa = flat * dtr
        phi = humr * ( t - 12.) + et
c
        a = sin(fa) * sin(dc)
        b = cos(fa) * cos(dc)
        	cosx = a + b * cos(phi)
        	if(abs(cosx).gt.1.) cosx=sign(1.,cosx)
        zenith = acos(cosx) / dtr
c
c calculate sunrise and sunset times --  at the ground...........
c see Explanatory Supplement to the Ephemeris (1961) pg 401......
c sunrise at height h metres is at...............................
		h=height*1000.
        chih = 90.83 + 0.0347 * sqrt(h)
c this includes corrections for horizontal refraction and........
c semi-diameter of the solar disk................................
        ch = cos(chih * dtr)
        cosphi = (ch -a ) / b
c if abs(secphi) > 1., sun does not rise/set.....................
c allow for sun never setting - high latitude summer.............
        secphi = 999999.
        if(cosphi.ne.0.) secphi = 1./cosphi
        sunset = 99.
        sunrse = 99.
        if(secphi.gt.-1.0.and.secphi.le.0.) return
c allow for sun never rising - high latitude winter..............
        sunset = -99.
        sunrse = -99.
        if(secphi.gt.0.0.and.secphi.lt.1.) return
c
        	cosx = cosphi
        	if(abs(cosx).gt.1.) cosx=sign(1.,cosx)
        phi = acos(cosx)
        et = et / humr
        phi = phi / humr
        sunrse = 12. - phi - et
        sunset = 12. + phi - et
        if(sunrse.lt.0.) sunrse = sunrse + 24.
        if(sunset.ge.24.) sunset = sunset - 24.
c
        return
        end
C
C
        SUBROUTINE MODA(IN,IYEAR,MONTH,IDAY,IDOY,NRDAYMO)
C-------------------------------------------------------------------
C CALCULATES DAY OF YEAR (IDOY, ddd) FROM YEAR (IYEAR, yy or yyyy), 
C MONTH (MONTH, mm) AND DAY OF MONTH (IDAY, dd) IF IN=0, OR MONTH 
C AND DAY FROM YEAR AND DAY OF YEAR IF IN=1. NRDAYMO is an output 
C parameter providing the number of days in the specific month.
C-------------------------------------------------------------------
        DIMENSION       MM(12)
        DATA            MM/31,28,31,30,31,30,31,31,30,31,30,31/

        IMO=0
        MOBE=0
c
c  leap year rule: years evenly divisible by 4 are leap years, except
c  years also evenly divisible by 100 are not leap years, except years 
c  also evenly divisible by 400 are leap years. The year 2000 therefore 
C  is a leap year. The 100 and 400 year exception rule
c     if((iyear/4*4.eq.iyear).and.(iyear/100*100.ne.iyear)) mm(2)=29
c  will become important again in the year 2100 which is not a leap 
C  year.
c
        mm(2)=28
        if(iyear/4*4.eq.iyear) mm(2)=29

        IF(IN.GT.0) GOTO 5
                mosum=0
                if(month.gt.1) then
                        do 1234 i=1,month-1 
1234                            mosum=mosum+mm(i)
                        endif
                idoy=mosum+iday
                nrdaymo=mm(month)
                RETURN

5       IMO=IMO+1
                IF(IMO.GT.12) GOTO 55
                MOOLD=MOBE
                nrdaymo=mm(imo)
                MOBE=MOBE+nrdaymo
                IF(MOBE.LT.IDOY) GOTO 5
55              MONTH=IMO
                IDAY=IDOY-MOOLD
        RETURN
        END             
C
C
           subroutine tcon(yr,mm,day,idn,rz,ig,rsn,nmonth)
c----------------------------------------------------------------
c input:        yr,mm,day       year(yyyy),month(mm),day(dd)
c               idn             day of year(ddd)
c output:       rz(3)           12-month-smoothed solar sunspot number
c               ig(3)           12-month-smoothed IG index
c               rsn             interpolation parameter
c               nmonth          previous or following month depending
c                               on day
c
c Requires I/O UNIT=12 to read the Rz12 and IG12 indices file IG_RZ.DAT 
c 
c rz(1) & ig(1) contain the indices for the month mm and rz(2) & ig(2)
c for the previous month (if day less than 15) or for the following
c month (otherwise). These indices are for the mid of the month. The
c indices for the given day are obtained by linear interpolation and
c are stored in rz(3) and ig(3).
c
c The indices file IG_RZ.DAT is structured as follows (values are 
c separated by comma): 
c   day, month, year of the last update of this file,
c   a blank line
c   start month, start year, end month, end year,
c   a blank line
c   the IG index for December of start year - 1 (needed for interpolation)
c   the 12 IG indices (13-months running mean) for start year, 
c   the 12 IG indices for the second year 
c       .. and so on until the end year,
c   the IG index for January of end year + 1 (needed for interpolation)
c   a blank line
c   the Rz index for December of start year - 1 (needed for interpolation)
c   the 12 Rz indices (13-months running mean) for the start year,
c   the 12 Rz indices for the second year 
c       .. and so on until the end year.
c   the Rz index for January of end year + 1 (needed for interpolation)
c 
c A negative Rz index means that the given index is the 13-months-
C running mean of the solar radio flux (F10.7). The close correlation 
C between (Rz)12 and (F10.7)12 is used to compute the (Rz)12 indices.
c
c An IG index of -111 indicates that no IG values are available for the
c time period. In this case a correlation function between (IG)12 and 
C (Rz)12 is used to obtain (IG)12.
c
c The computation of the 13-month-running mean for month M requires the
c indices for the six months preceeding M and the six months following 
C M (month: M-6, ..., M+6). To calculate the current running mean one 
C therefore requires predictions of the indix for the next six months. 
C Starting from six months before the UPDATE DATE (listed at the top of 
c the file) and onward the indices are therefore based on indices 
c predictions.
c----------------------------------------------------------------

           integer	yr, mm, day, iflag, iyst, iyend,iymst
           integer	imst,iymend
           real		ionoindx(806),indrz(806)
           real		ig(3),rz(3)
           logical	mess
           
           common 	/iounit/konsol,mess

           save     ionoindx,indrz,iflag,iyst,iymst,iymend,imst

        if(iflag.eq.0) then      
            open(unit=12,file='ig_rz.dat',status='old')

c-web- special for web version
c            open(unit=12,file=
c     *         '/usr/local/etc/httpd/cgi-bin/models/IRI/ig_rz.dat',
c     *         status='old')

c Read the update date, the start date and the end date (mm,yyyy), and
c get number of data points to read.

            read(12,*) iupd,iupm,iupy
            read(12,*) imst,iyst,imend,iyend
            iymst=iyst*100+imst
            iymend=iyend*100+imend

c inum_vals= 12-imst+1+(iyend-iyst-1)*12 +imend + 2
c 1st year \ full years       \last y\ before & after

            inum_vals= 3-imst+(iyend-iyst)*12 +imend

c read all the IG12 (ionoindx) and Rz12 (indrz) values

            read(12,*) (ionoindx(i),i=1,inum_vals)
            read(12,*) (indrz(i),i=1,inum_vals)
            do 1 jj=1,inum_vals
                rrr=indrz(jj)
                if(rrr.lt.0.0) then
                    covr=abs(rrr)
                    rrr=33.52*sqrt(covr+85.12)-408.99
                    if(rrr.lt.0.0) rrr=0.0
                    indrz(jj)=rrr
                    endif
                if(ionoindx(jj).gt.-90.) goto 1
                    zi=-12.349154+(1.4683266-2.67690893e-03*rrr)*rrr
                    if(zi.gt.274.0) zi=274.0
                    ionoindx(jj)=zi
1               continue
            close(unit=12)
            iflag = 1
        endif

        iytmp=yr*100+mm
        if (iytmp .lt. iymst .or. iytmp .gt. iymend) then
               if(mess) write(konsol,8000) iytmp,iymst,
     &                                            iymend
 8000          format(1x,I10,'** OUT OF RANGE **'/,5x,
     &  'The file IG_RZ.DAT which contains the indices Rz12',
     &  ' and IG12'/5x,'currently only covers the time period',
     &  ' (yymm) : ',I6,'-',I6)
               nmonth=-1
               return
               endif

c       num=12-imst+1+(yr-iyst-1)*12+mm+1
        num=2-imst+(yr-iyst)*12+mm

        rz(1)=indrz(num)
        ig(1)=ionoindx(num)
        midm=15
        if(mm.eq.2) midm=14
        call MODA(0,yr,mm,midm,idd1,nrdaym)

        if(day.lt.midm) goto 1926
c day is at or after mid of month
                imm2=mm+1
                if(imm2.gt.12) then
                        imm2=1
                        iyy2=yr+1
                        idd2=380            ! =365+15 mid-January
c               if((yr/4*4.eq.yr).and.(yr/100*100.ne.yr)) idd2=381
                        if(yr/4*4.eq.yr) idd2=381
                else
                        iyy2=yr
                        midm=15
                        if(imm2.eq.2) midm=14
                        call MODA(0,iyy2,imm2,midm,IDD2,nrdaym)
                endif
                rz(2)=indrz(num+1)
                ig(2)=ionoindx(num+1)
                rsn=(idn-idd1)*1./(idd2-idd1)                
                rz(3)=rz(1)+(rz(2)-rz(1))*rsn
                ig(3)=ig(1)+(ig(2)-ig(1))*rsn
                goto 1927
1926            imm2=mm-1
                if(imm2.lt.1) then
                        imm2=12
                        idd2=-16
                        iyy2=yr-1
                else
                        iyy2=yr
                        midm=15
                        if(imm2.eq.2) midm=14
                        call MODA(0,iyy2,imm2,midm,IDD2,nrdaym)
                endif
                rz(2)=indrz(num-1)
                ig(2)=ionoindx(num-1)
                rsn=(idn-idd2)*1./(idd1-idd2)
                rz(3)=rz(2)+(rz(1)-rz(2))*rsn
                ig(3)=ig(2)+(ig(1)-ig(2))*rsn

1927    nmonth=imm2
            return
            end
C
C
        SUBROUTINE APF_ONLY(IYYYY,IMN,ID,F107D,F107PD,F107_81,F107_365,
     *        IAPDA)
c-----------------------------------------------------------------------
c Finds daily F10.7, daily Ap, and 81-day and 365-day F10.7 index: 
c
c    INPUTS: 	IYYYY (yyyy)	year 
c 				IMN (mm)		month 
c				ID (dd)			day 
c    OUTPUT:    F107D			F10.7 index for the day (adjusted 
c								to 1AU)
C               F107PD  		F10.7 index for one day prior (used in MSIS)
c				F107_81			F10.7 average over 3 solar rotations
c                               (81 days, centered on the current day) 
c               F107_365        F10.7 12-month running mean
c				IAPDA			Daily Ap
c 
c Using APF107.DAT file (see subroutine APF) on UNIT=13. 
c
c Is used for vdrift and foeedi.
c
c If date is outside the range of indices file than F107D=F107_81=-11.1  
c-----------------------------------------------------------------------

        DIMENSION iiap(8),lm(12)
        LOGICAL mess

        common /iounit/konsol,mess

        DATA LM/31,28,31,30,31,30,31,31,30,31,30,31/

        IYBEG=1958
        if(iyyyy.lt.IYBEG) goto 21   ! APF107.DAT starts at Jan 1, 1958

        Open(13,FILe='apf107.dat',
c-web-sepcial vfor web version
C      OPEN(13,FILE='/usr/local/etc/httpd/cgi-bin/models/IRI/apf107.dat',
     *    ACCESS='DIRECT',RECL=55,FORM='FORMATTED',STATUS='OLD')

        is=0
        do i=IYBEG,iyyyy-1
            nyd=365
            if(i/4*4.eq.i) nyd=366	! leap year
            IS=IS+nyd
            enddo

        lm(2)=28
        if(iyyyy/4*4.eq.iyyyy) lm(2)=29	  ! leap year
                do i=1,IMN-1
              IS=IS+LM(i)
              ENDDO
        
        IS=IS+ID

        READ(13,10,REC=IS,ERR=21) JY,JMN,JD,iiap,iapda,IR,F107D,F107_81,
     *        F365
 		
 		if(F107_81.lt.-4.) F107_81=F107D
 		if(F107_365.lt.-4.) F107_365=F107D

        F107PD=F107D
		if(IS.gt.1) READ(13,10,REC=IS-1,ERR=21) JY,JMN,JD,iiap,
     *  	idum1,idum2,F107PD,fdum1,fdum2
                 
10      FORMAT(3I3,9I3,I3,3F5.1)
        CLOSE(13)
        goto 20

21      if(mess) write(konsol,100)
100     format(1X,'APF_ONLY: Date is outside range of F10.7D indices',
     &    ' file (F10.7D = F10.7_81 = F10.7RM12).')
        F107D = -11.1
        F107_81 = -11.1
        F107_365 = -11.1
     
20    RETURN
      END

C
C
        SUBROUTINE APF(IYYYY,IMN,ID,HOUR,IAP)
c-----------------------------------------------------------------------
c Finds 3-hourly Ap indices for IRI-STORM model
c    INPUTS: 	IYYYY 			year (yyyy)
c 				IMN 			month (mm)
c				ID 				day (dd)	
c				HOUR			UT in decimal hours
c    OUTPUT:    IAP(1:13)		3-hourly Ap index
c								IAP(13) Ap index for current UT
c								IAP(1) AP index for UT-39 hours.
c
c Reads APF107.DAT file (on UNIT=13) that is structured as follows:
c 		JY(I3),JMN(I3),JD(I3)	year, month, day 
c		IIAP(8)	(8I3)			3-hour Ap indices for the UT intervals 
c								(0-3),)3-6),)6-9), .., )18-21),)21-24(
c		IAPD (I3)				daily Ap
c		IR (I3)					sunspot number for the day (empty)
c		F107 (F5.1)				F10.7 radio flux for the day
c		F107_81 (F5.1)			81-day average of F10.7 radio flux 
c       F107_365 (F5.1)         365-day average of F10.7 centered on 
c                               the date of interest. At start and end  
c								of index file it takes all available  
c                               indices, e.g. for the first date the 
c                               average is only over 40 F10.7 values  
c                               and over 41 values on the 2nd date.  
c
c If date is outside the range of the Ap indices file than IAP(1)=-5  
c-----------------------------------------------------------------------

        DIMENSION iiap(8),iap(13),lm(12)
        LOGICAL	mess
        COMMON /iounit/konsol,mess
        DATA LM/31,28,31,30,31,30,31,31,30,31,30,31/

        IYBEG=1958
       
        do i=1,8
              iap(i)=-1
              enddo

        if(iyyyy.lt.IYBEG) goto 21   ! file starts at Jan 1, 1958

        OPEN(13,FILe='apf107.dat',
c-web-sepcial vfor web version
c      OPEN(13,FILE='/usr/local/etc/httpd/cgi-bin/models/IRI/apf107.dat',
     *    ACCESS='DIRECT',RECL=55,FORM='FORMATTED',STATUS='OLD')
                
        is=0
        if(iyyyy.gt.IYBEG) then
            do i=IYBEG,iyyyy-1
                nyd=365
                if(i/4*4.eq.i) nyd=366
                IS=IS+nyd
                enddo
            endif   

        lm(2)=28
        if(iyyyy/4*4.eq.iyyyy) lm(2)=29
        do i=1,IMN-1
              IS=IS+LM(i)
              ENDDO

        IS=IS+ID

        ihour=int(hour/3.)+1
        if(ihour.gt.8) ihour=8

        if(is*8+ihour.lt.13) goto 21   ! at least 13 indices available	
        READ(13,10,REC=IS,ERR=21) JY,JMN,JD,iiap,iapd,IR,F,F81,F365
        do i9=1,8
        	if(iiap(i9).lt.-2) goto 21
        	enddo
        j1=13-ihour
        do i=1,ihour
           iap(j1+i)=iiap(i)
           enddo
        iss=is-1
        READ(13,10,REC=ISS,ERR=21) JY,JMN,JD,iiap,iapd,IR,F,F81,F365
        do i9=1,8
        	if(iiap(i9).lt.-2) goto 21
        	enddo
        if(ihour.gt.4) then
              do i=1,j1
                iap(i)=iiap(8-j1+i)
                enddo
        else           
             j2=5-ihour
             do i=1,8
                iap(j2+i)=iiap(i)
                enddo
             iss=is-2
             READ(13,10,REC=ISS,ERR=21) JY,JMN,JD,iiap,iapd,IR,F,F81,F365
        	 do i9=1,8
        		if(iiap(i9).lt.-2) goto 21
        		enddo
             do i=1,j2
                iap(i)=iiap(8-j2+i)
                enddo
        endif         
  10    FORMAT(3I3,9I3,I3,3F5.1)
        CLOSE(13)
        goto 20
        
21      if(mess) write(konsol,100)
100     format(1X,'Date is outside range of Ap indices file.',
     &     ' STORM model is turned off.')
        IAP(1)=-5
      
20    RETURN
      END
C
C
C----------------------STORM MODEL --------------------------------
C
      SUBROUTINE CONVER(rga,rgo,rgma)

C     This subroutine converts a geographic latitude and longitude
C     location to a corrected geomagnetic latitude.
C
C     INPUT: 
C       geographic latitude   -90. to +90.
C       geographic longitude  0. to 360. positive east from Greenwich.
C
C     OUTPUT:
C       corrected geomagnetic latitude	-90. to +90.


      DIMENSION CORMAG(20,91)      
      DATA ((CORMAG(i,j),i=1,20),j=1,31)/
     +163.68,163.68,163.68,163.68,163.68,163.68,
     +163.68,163.68,163.68,163.68,163.68,163.68,163.68,163.68,
     +163.68,163.68,163.68,163.68,163.68,163.68,162.60,163.12,
     +163.64,164.18,164.54,164.90,165.16,165.66,166.00,165.86,
     +165.20,164.38,163.66,162.94,162.42,162.00,161.70,161.70,
     +161.80,162.14,161.20,162.18,163.26,164.44,165.62,166.60,
     +167.42,167.80,167.38,166.82,166.00,164.66,163.26,162.16,
     +161.18,160.40,159.94,159.80,159.98,160.44,159.80,161.14,
     +162.70,164.50,166.26,167.90,169.18,169.72,169.36,168.24,
     +166.70,164.80,162.90,161.18,159.74,158.60,157.94,157.80,
     +157.98,158.72,158.40,160.10,162.02,164.28,166.64,169.00,
     +170.80,171.72,171.06,169.46,167.10,164.64,162.18,160.02,
     +158.20,156.80,156.04,155.80,156.16,157.02,157.00,158.96,
     +161.24,163.86,166.72,169.80,172.42,173.72,172.82,170.34,
     +167.30,164.22,161.34,158.74,156.60,155.00,154.08,153.90,
     +154.36,155.36,155.50,157.72,160.36,163.32,166.60,170.20,
     +173.70,175.64,174.18,170.80,167.10,163.56,160.24,157.36,
     +154.96,153.10,152.08,151.92,152.46,153.76,154.10,156.52,
     +159.36,162.52,166.24,170.30,174.62,177.48,175.04,170.82,
     +166.60,162.70,159.02,155.88,153.22,151.20,150.08,149.92,
     +150.64,152.20,152.80,155.32,158.28,161.70,165.58,170.00,
     +174.84,178.46,175.18,170.38,165.80,161.64,157.80,154.38,
     +151.52,149.30,148.18,148.02,148.92,150.60,151.40,154.08,
     +157.18,160.68,164.78,169.40,174.34,177.44,174.28,169.44,
     +164.70,160.34,156.30,152.78,149.72,147.40,146.18,146.04,
     +147.12,149.04,150.10,152.88,156.00,159.58,163.78,168.50,
     +173.28,175.60,172.86,168.14,163.40,158.98,154.88,151.10,
     +147.98,145.50,144.18,144.14,145.40,147.48,148.80,151.68,
     +154.88,158.48,162.68,167.40,171.76,173.60,171.12,166.68,
     +162.00,157.48,153.28,149.50,146.18,143.50,142.18,142.24,
     +143.68,145.98,147.50,150.54,153.68,157.28,161.42,166.10,
     +170.10,171.48,169.22,164.98,160.40,155.88,151.68,147.80,
     +144.34,141.60,140.18,140.26,141.98,144.62,146.30,149.34,
     +152.48,155.98,160.08,164.60,168.34,169.38,167.20,163.18,
     +158.60,154.18,149.98,146.02,142.54,139.70,138.18,138.46,
     +140.26,143.16,145.10,148.14,151.18,154.60,158.68,163.10,
     +166.48,167.28,165.18,161.32,156.90,152.48,148.28,144.32,
     +140.74,137.80,136.22,136.48,138.64,141.76,143.90,146.98,
     +149.98,153.30,157.24,161.40,164.52,165.16,162.86,159.42,
     +155.00,150.68,146.48,142.52,138.94,135.90,134.22,134.68,
     +137.02,140.40,142.70,145.84,148.76,151.92,155.74,159.70,
     +162.52,162.96,160.98,157.42,153.10,148.84,144.68,140.82,
     +137.20,134.00,132.32,132.80,135.42,139.10,141.60,144.74,
     +147.46,150.52,154.20,158.00,160.46,160.76,158.86,155.36,
     +151.20,146.94,142.88,139.02,135.40,132.10,130.32,131.00,
     +133.80,137.74,140.50,143.58,146.24,149.12,152.60,156.20,
     +158.40,158.66,156.76,153.36,149.30,145.04,141.08,137.30,
     +133.60,130.30,128.42,129.12,132.28,136.44,139.30,142.48,
     +144.94,147.64,150.48,154.30,156.34,156.36,154.56,151.26,
     +147.30,143.14,139.20,135.50,131.90,128.40,126.52,127.32,
     +130.76,135.18,138.20,141.28,143.72,146.24,149.26,152.40,
     +154.24,154.16,152.36,149.16,145.30,141.24,137.30,133.70,
     +130.10,126.60,124.62,125.54,129.16,133.92,137.10,140.18,
     +142.42,144.66,147.62,150.50,152.18,151.96,150.16,147.10,
     +143.30,139.24,135.50,131.90,128.36,124.80,122.72,123.74,
     +127.64,132.62,135.90,139.02,141.12,143.18,145.92,148.60,
     +149.98,149.76,148.04,145.00,141.20,137.30,133.60,130.10,
     +126.60,123.00,120.86,121.96,126.12,131.36,134.80,137.88,
     +139.80,141.68,144.08,146.60,147.88,147.56,145.84,142.90,
     +139.20,135.30,131.70,128.28,124.86,121.30,118.96,120.18,
     +124.70,130.16,133.60,136.72,138.48,140.10,142.38,144.60,
     +145.72,145.34,143.64,140.80,137.10,133.30,129.72,126.48,
     +123.10,119.50,117.16,118.48,123.18,128.86,132.40,135.42,
     +137.08,138.50,140.54,142.60,143.52,143.06,141.44,138.70,
     +135.10,131.30,127.82,124.58,121.40,117.70,115.26,116.70,
     +121.66,127.60,131.20,134.22,135.66,136.82,138.70,140.60,
     +141.36,140.86,139.24,136.50,133.00,129.30,125.92,122.78,
     +119.60,116.00,113.40,114.92,120.16,126.30,130.00,132.92,
     +134.24,135.14,136.80,138.60,139.16,138.64,137.12,134.40,
     +130.90,127.20,123.92,120.96,117.90,114.20,111.56,113.12,
     +118.64,124.90,128.70,131.56,132.74,133.44,134.90,136.50,
     +137.00,136.36,134.82,132.30,128.70,125.16,121.94,119.06,
     +116.10,112.50,109.70,111.42,117.14,123.60,127.30,130.16,
     +131.22,131.66,133.00,134.50,134.80,134.14,132.62,130.14,
     +126.60,123.06,119.94,117.16,114.30,110.70,107.80,109.64,
     +115.62,122.24,125.90,128.76,129.62,129.96,131.06,132.40,
     +132.60,131.86,130.42,128.00,124.50,120.96,117.96,115.26,
     +112.54,108.90,105.94,107.86,114.02,120.84/

      DATA ((CORMAG(i,j),i=1,20),j=32,61)/
     +124.05,126.79,
     +127.55,127.83,128.90,130.21,130.41,129.71,128.33,125.96,
     +122.49,118.96,115.97,113.26,110.52,106.89,104.01,106.00,
     +112.21,119.06,122.19,124.82,125.48,125.69,126.73,128.03,
     +128.22,127.55,126.23,123.92,120.47,116.97,113.97,111.26,
     +108.50,104.89,102.08,104.14,110.41,117.29,120.34,122.85,
     +123.40,123.56,124.57,125.84,126.03,125.40,124.14,121.88,
     +118.46,114.97,111.98,109.26,106.48,102.88,100.15,102.28,
     +108.60,115.51,118.49,120.88,121.33,121.42,122.40,123.65,
     +123.84,123.24,122.04,119.83,116.45,112.97,109.98,107.26,
     +104.46,100.87,098.22,100.42,106.79,113.74,116.63,118.91,
     +119.26,119.29,120.24,121.47,121.65,121.09,119.95,117.79,
     +114.43,110.98,107.99,105.26,102.44,098.87,096.29,098.56,
     +104.98,111.96,114.78,116.94,117.19,117.15,118.07,119.28,
     +119.46,118.93,117.86,115.75,112.42,108.98,106.00,103.26,
     +100.42,096.86,094.36,096.70,103.18,110.19,112.93,114.97,
     +115.12,115.02,115.91,117.09,117.27,116.78,115.76,113.71,
     +110.41,106.98,104.00,101.26,098.40,094.85,092.43,094.84,
     +101.37,108.41,111.07,113.00,113.04,112.88,113.74,114.91,
     +115.08,114.62,113.67,111.67,108.39,104.99,102.01,099.26,
     +096.38,092.85,090.51,092.97,099.56,106.64,109.22,111.03,
     +110.97,110.75,111.58,112.72,112.89,112.47,111.57,109.63,
     +106.38,102.99,100.01,097.26,094.36,090.84,088.58,091.11,
     +097.75,104.86,107.37,109.06,108.90,108.61,109.41,110.53,
     +110.70,110.31,109.48,107.59,104.37,100.99,098.02,095.26,
     +092.34,088.83,086.65,089.25,095.95,103.09,105.51,107.09,
     +106.83,106.48,107.25,108.35,108.51,108.16,107.39,105.55,
     +102.35,099.00,096.03,093.26,090.32,086.83,084.72,087.39,
     +094.14,101.31,103.66,105.12,104.76,104.34,105.08,106.16,
     +106.32,106.00,105.29,103.50,100.34,097.00,094.03,091.26,
     +088.30,084.82,082.79,085.53,092.33,099.54,101.81,103.15,
     +102.68,102.21,102.92,103.97,104.13,103.85,103.20,101.46,
     +098.33,095.00,092.04,089.26,086.28,082.81,080.86,083.67,
     +090.52,097.76,099.95,101.18,100.61,100.07,100.75,101.79,
     +101.94,101.69,101.10,099.42,096.31,093.01,090.04,087.26,
     +084.26,080.81,078.93,081.81,088.72,095.99,098.10,099.21,
     +098.54,097.94,098.59,099.60,099.75,099.54,099.01,097.38,
     +094.30,091.01,088.05,085.26,082.24,078.80,077.00,079.95,
     +086.91,094.21,096.25,097.24,096.47,095.81,096.43,097.41,
     +097.56,097.39,096.92,095.34,092.29,089.01,086.06,083.26,
     +080.22,076.79,075.07,078.09,085.10,092.43,094.39,095.27,
     +094.40,093.67,094.26,095.23,095.37,095.23,094.82,093.30,
     +090.27,087.02,084.06,081.26,078.20,074.79,073.14,076.23,
     +083.30,090.66,092.54,093.30,092.32,091.54,092.10,093.04,
     +093.18,093.08,092.73,091.26,088.26,085.02,082.07,079.26,
     +076.18,072.78,071.21,074.37,081.49,088.88,090.69,091.33,
     +090.25,089.40,089.93,090.85,090.99,090.92,090.63,089.21,
     +086.25,083.02,080.07,077.26,074.16,070.77,069.28,072.51,
     +079.68,087.11,088.83,089.36,088.18,087.27,087.77,088.67,
     +088.80,088.77,088.54,087.17,084.23,081.03,078.08,075.26,
     +072.14,068.77,067.35,070.65,077.87,085.33,086.98,087.39,
     +086.11,085.13,085.60,086.48,086.61,086.61,086.45,085.13,
     +082.22,079.03,076.09,073.26,070.12,066.76,065.42,068.79,
     +076.07,083.56,085.13,085.42,084.04,083.00,083.44,084.29,
     +084.42,084.46,084.35,083.09,080.21,077.03,074.09,071.26,
     +068.10,064.75,063.49,066.93,074.26,081.78,083.27,083.45,
     +081.96,080.86,081.27,082.11,082.23,082.30,082.26,081.05,
     +078.19,075.04,072.10,069.26,066.08,062.75,061.57,065.06,
     +072.45,080.01,081.42,081.48,079.89,078.73,079.11,079.92,
     +080.04,080.15,080.16,079.01,076.18,073.04,070.10,067.26,
     +064.06,060.74,059.64,063.20,070.64,078.23,079.57,079.51,
     +077.82,076.59,076.94,077.73,077.85,077.99,078.07,076.97,
     +074.17,071.04,068.11,065.26,062.04,058.73,057.71,061.34,
     +068.84,076.46,077.71,077.54,075.75,074.46,074.78,075.55,
     +075.66,075.84,075.98,074.93,072.15,069.05,066.12,063.26,
     +060.02,056.73,055.78,059.48,067.03,074.68,075.86,075.57,
     +073.68,072.32,072.61,073.36,073.47,073.68,073.88,072.88,
     +070.14,067.05,064.12,061.26,058.00,054.72,053.85,057.62,
     +065.22,072.91,074.01,073.60,071.60,070.19,070.45,071.17,
     +071.28,071.53,071.79,070.84,068.13,065.05,062.13,059.26,
     +055.98,052.71,051.92,055.76,063.41,071.13,072.15,071.63,
     +069.53,068.05,068.28,068.99,069.09,069.37,069.69,068.80,
     +066.11,063.06,060.13,057.26,053.96,050.71,049.99,053.90,
     +061.61,069.36,070.30,069.66,067.46,065.92,066.12,066.80,
     +066.90,067.22,067.60,066.76,064.10,061.06,058.14,055.26,
     +051.94,048.70,048.06,052.04,059.80,067.58/

      DATA ((CORMAG(i,j),i=1,20),j=62,91)/
     +067.70,067.06,
     +065.08,063.72,063.98,064.60,064.80,065.12,065.60,064.86,
     +062.40,059.26,056.24,053.18,049.84,046.60,046.12,050.12,
     +057.52,064.80,064.90,064.42,062.70,061.62,061.78,062.40,
     +062.60,063.04,063.58,063.00,060.60,057.46,054.42,051.18,
     +047.70,044.60,044.22,048.02,055.06,061.92,062.10,061.72,
     +060.32,059.50,059.68,060.20,060.46,060.94,061.58,061.00,
     +058.70,055.66,052.52,049.18,045.60,042.50,042.22,046.00,
     +052.60,058.98,059.20,059.18,058.12,057.32,057.48,058.00,
     +058.30,058.84,059.48,059.04,056.90,053.86,050.62,047.10,
     +043.50,040.50,040.28,043.98,050.22,056.18,056.40,056.64,
     +055.84,055.20,055.38,055.80,056.16,056.84,057.48,057.04,
     +055.10,052.06,048.70,045.10,041.40,038.40,038.28,041.88,
     +047.94,053.44,053.70,054.14,053.56,053.10,053.24,053.70,
     +054.06,054.74,055.38,055.14,053.20,050.26,046.80,043.10,
     +039.34,036.40,036.38,039.96,045.56,050.84,051.10,051.70,
     +051.36,051.00,051.14,051.50,051.96,052.64,053.38,053.08,
     +051.30,048.36,044.90,041.02,037.24,034.40,034.38,037.86,
     +043.28,048.20,048.50,049.26,049.18,048.90,049.04,049.40,
     +049.86,050.64,051.28,051.08,049.40,046.46,042.98,039.02,
     +035.14,032.40,032.48,035.72,041.00,045.70,046.00,046.96,
     +046.98,046.80,046.94,047.30,047.76,048.54,049.28,049.08,
     +047.40,044.56,041.08,037.02,033.14,030.40,030.58,033.84,
     +038.72,043.20,043.50,044.62,044.80,044.80,044.94,045.20,
     +045.76,046.54,047.18,046.98,045.50,042.66,039.08,035.02,
     +031.14,028.40,028.58,031.82,036.52,040.80,041.20,042.32,
     +042.54,042.70,042.84,043.20,043.66,044.44,045.08,044.98,
     +043.50,040.76,037.08,033.04,029.04,026.40,026.68,029.82,
     +034.34,038.40,038.80,040.12,040.60,040.70,040.84,041.10,
     +041.62,042.34,042.98,042.88,041.50,038.76,035.18,031.04,
     +027.14,024.50,024.78,027.70,032.14,036.06,036.50,037.88,
     +038.50,038.68,038.84,039.10,039.56,040.34,040.88,040.82,
     +039.40,036.76,033.18,029.12,025.14,022.50,022.88,025.90,
     +029.96,033.86,034.30,035.68,036.42,036.68,036.84,037.10,
     +037.56,038.24,038.88,038.72,037.40,034.76,031.18,027.12,
     +023.14,020.60,020.98,023.90,027.88,031.66,032.10,033.58,
     +034.32,034.68,034.84,035.10,035.56,036.24,036.78,036.62,
     +035.30,032.72,029.18,025.14,021.24,018.70,019.08,021.90,
     +025.88,029.42,029.90,031.48,032.32,032.68,032.84,033.10,
     +033.56,034.22,034.68,034.42,033.20,030.72,027.28,023.22,
     +019.34,016.80,017.24,020.00,023.78,027.32,027.70,029.38,
     +030.24,030.68,030.94,031.20,031.66,032.22,032.58,032.32,
     +031.10,028.62,025.28,021.32,017.48,015.00,015.38,018.18,
     +021.80,025.22,025.70,027.28,028.24,028.78,029.04,029.30,
     +029.66,030.22,030.50,030.22,029.00,026.62,023.30,019.42,
     +015.64,013.10,013.54,016.28,019.80,023.12,023.60,025.24,
     +026.24,026.78,027.14,027.40,027.76,028.22,028.40,028.12,
     +026.80,024.52,021.30,017.52,013.78,011.30,011.74,014.48,
     +017.90,021.12,021.60,023.24,024.34,024.88,025.24,025.50,
     +025.86,026.22,026.40,025.98,024.70,022.48,019.40,015.72,
     +012.04,009.50,009.94,012.58,016.02,019.12,019.60,021.24,
     +022.34,022.98,023.34,023.70,024.00,024.30,024.40,023.88,
     +022.60,020.48,017.52,014.00,010.34,007.80,008.18,010.88,
     +014.22,017.18,017.60,019.34,020.44,021.16,021.54,021.90,
     +022.16,022.40,022.32,021.78,020.60,018.48,015.62,012.20,
     +008.68,006.00,006.44,009.18,012.42,015.28,015.80,017.44,
     +018.54,019.26,019.74,020.10,020.30,020.50,020.32,019.72,
     +018.50,016.54,013.84,010.68,007.14,004.40,004.74,007.58,
     +010.74,013.48,014.00,015.54,016.74,017.46,017.94,018.30,
     +018.50,018.58,018.32,017.72,016.50,014.64,012.24,009.18,
     +005.84,002.90,003.30,006.16,009.14,011.84,012.30,013.78,
     +014.94,015.66,016.24,016.50,016.70,016.70,016.42,005.78,
     +014.60,012.90,010.66,007.86,004.88,001.60,001.72,004.96,
     +007.84,010.24,010.70,012.14,013.24,013.96,014.44,014.80,
     +014.90,014.88,014.52,013.92,012.80,011.30,009.28,006.94,
     +004.32,001.80,001.94,004.34,006.78,008.94,009.40,010.58,
     +011.64,012.36,012.74,013.10,013.20,013.08,012.72,012.12,
     +011.10,009.86,008.30,006.50,004.60,003.10,003.16,004.50,
     +006.20,007.90,008.40,009.42,010.14,010.76,011.14,011.40,
     +011.40,011.38,011.02,010.46,009.70,008.72,007.64,006.46,
     +005.42,004.60,004.70,005.34,006.24,007.36,007.90,008.46,
     +008.92,009.28,009.54,009.70,009.70,009.68,009.42,009.06,
     +008.60,008.08,007.56,007.02,006.56,006.30,006.30,006.52,
     +006.96,007.38,008.15,008.15,008.15,008.15,008.15,008.15,
     +008.15,008.15,008.15,008.15,008.15,008.15,008.15,008.15,
     +008.15,008.15,008.15,008.15,008.15,008.15/

C     Data Input      
      rlan = rga
      rlo = rgo      
      
C     From "normal" geographic latitude 
C     to angle from South Pole.       
      rla = rlan + 90

      IF (rlo .EQ. 360) THEN
      	rlo = 0
        END IF

C     PROXIMITY

C     coefficients of the latitudinal points		
      LA1 = (INT(rla/2)+1)
      LA2 = LA1 + 1
      if(la2.gt.91) la2=91

C     coefficients of the longitudinal points		
      LO1 = (INT(rlo/18)+1)
corr      LO2 = LO1 + 1
      LO2 = MOD(LO1,20) + 1 

C     Four points of Geomagnetic Coordinates
      gm1 = CORMAG(LO1,LA1)
      gm2 = CORMAG(LO1,LA2) 
      gm3 = CORMAG(LO2,LA1)
      gm4 = CORMAG(LO2,LA2)

C     latitudinal points		
C      X1 = ABS(rla - (INT(rla)))                        
C      X2 = 2. - X1
	  x = (rla/2.0 - (INT(rla/2.0)))

C     longitudinal points		
C      Y1 = ABS(rlo - (INT(rlo)))                        
C      Y2 = 18. - Y1
      y =(rlo/18.0 - (INT(rlo/18.0))) 

C     X AND Y VALUES
C      x = X1 / (X1 + X2)
C      y = Y1 / (Y1 + Y2)

C     INTERPOLATION
      gmla = gm1 * (1 - x) * (1 - y) + gm2 * (1 - y) * (x) + gm3 * (y)
     1 * (1 - x) + gm4 * (x) * (y)

C     OUTPUT OF THE PROGRAM
C     From corrected geomagnetic latitude from North Pole
C     to "normal"  geomagnetic latitude.       
      rgma = 90. - gmla

      END
c
c
      SUBROUTINE STORM(ap,rga,rgo,coor,rgma,ut,doy,cf)
C----------------------------------------------------------------------
C      Fortran code to obtain the foF2 storm-time correction factor at 
C      a given location and time, using the current and the 12 previous
C      ap values as input.
C
C      ap ---> (13 elements integer array). Array with the preceeding
C              13 value of the 3-hourly ap index. The 13th value
C              in the array will contain the ap at the UT of interest,
C              the 12th value will contain the 1st three hourly interval
C              preceeding the time of interest, and so on to the 1st
C              ap value at the earliest time.
C     coor --> (integer). If coor = 2, rga should contain the 
C                         geomagnetic latitude.
C                         If coor = 1, rga should contain the 
C                         geographic latitude.
C     rga ---> (real, -90 to 90) geographic or geomagnetic latitude.
C     rgo ---> (real, 0 to 360, positive east from Greenwich.)
C                           geographic longitude, only used if coor=1.
C     ut  ---> (integer, hours 00 to 23) Universal Time of interest.
C     doy ---> (integer, 1 to 366)Day of the year.
C     cf  ---> (real) The output; the storm-time correction factor used
C              to scale foF2, foF2 * cf.
C     rgma --> corrected magnetic latitude calculated from rga and rgo
C
C     This model and computer code was developed by E. Araujo-Pradere,
C     T. Fuller-Rowell and M. Condrescu, SEC, NOAA, Boulder, USA
C     Ref: 
C     T. Fuller-Rowell, E. Araujo-Pradere, and M. Condrescu, An 
C       Empirical Ionospheric Storm-Time Ionospheric Correction Model,
C       Adv. Space Res. 8, 8, 15-24, 2000.
C----------------------------------------------------------------------
C     DIMENSIONS AND COEFFICIENTS VALUES

      DIMENSION c4(20)
      DATA c4/0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,
     +0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,
     +0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00/

      DIMENSION c3(20)
      DATA c3/0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,-9.44E-12,
     +0.00E+00,3.04E-12,0.00E+00,9.32E-12,-1.07E-11,0.00E+00,0.00E+00,
     +0.00E+00,1.09E-11,0.00E+00,0.00E+00,0.00E+00,0.00E+00,-1.01E-11/

      DIMENSION c2(20)
      DATA c2/1.16E-08,0.00E+00,0.00E+00,-1.46E-08,0.00E+00,9.86E-08,
     +2.25E-08,-1.67E-08,-1.62E-08,-9.42E-08,1.17E-07,4.32E-08,3.97E-08,
     +3.13E-08,-8.04E-08,3.91E-08,2.58E-08,3.45E-08,4.76E-08,1.13E-07/

      DIMENSION c1(20)
      DATA c1/-9.17E-05,-1.37E-05,0.00E+00,7.14E-05,0.00E+00,-3.21E-04,
     +-1.66E-04,-4.10E-05,1.36E-04,2.29E-04,-3.89E-04,-3.08E-04,
     +-2.81E-04,-1.90E-04,4.76E-05,-2.80E-04,-2.07E-04,-2.91E-04,
     +-3.30E-04,-4.04E-04/

      DIMENSION c0(20)
      DATA c0/1.0136E+00,1.0478E+00,1.00E+00,1.0258E+00,1.00E+00,
     +1.077E+00,1.0543E+00,1.0103E+00,9.9927E-01,9.6876E-01,1.0971E+00,
     +1.0971E+00,1.0777E+00,1.1134E+00,1.0237E+00,1.0703E+00,1.0248E+00,
     +1.0945E+00,1.1622E+00,1.1393E+00/

      DIMENSION fap(36)
      DATA fap/0.,0.,0.037037037,0.074074074,0.111111111,0.148148148,
     10.185185185,0.222222222,0.259259259,0.296296296,0.333333333,
     20.37037037,0.407407407,0.444444444,0.481481481,0.518518519,
     30.555555556,0.592592593,0.62962963,0.666666667,0.703703704,
     40.740740741,0.777777778,0.814814815,0.851851852,0.888888889,
     50.925925926,0.962962963,1.,0.66666667,0.33333334,0.,0.333333,
     60.666666,1.,0.7/

      integer code(8,6)
      data code/3,4,5,4,3,2,1,2,3,2,1,2,3,4,5,4,8,7,6,7,8,9,10,9,
     *13,12,11,12,13,14,15,14,18,17,16,17,18,19,20,19,18,17,16,17,
     *18,19,20,19/

      INTEGER ape(39)
      INTEGER ap(13)
      INTEGER ut,doy,dayno,coor,s1,s2,l1,l2
      REAL rgma, rap, rga, rgo, rs, rl

C      CALLING THE PROGRAM TO CONVERT TO GEOMAGNETIC COORDINATES

       IF (coor .EQ. 1) THEN

           CALL CONVER (rga,rgo,rgma)

       ELSE IF (coor .EQ. 2) THEN
                rgma = rga

       ELSE

          WRITE (6,*)' '
          WRITE (6,*)' '
          WRITE (6,*)'   Wrong Coordinates Selection -------- >>', coor
          WRITE (6,*)' '
          GOTO 100
       ENDIF

C FROM 3-HOURLY TO HOURLY ap (New, interpolates between the three hourly ap values)

       ape(1)=ap(1)
       ape(2)=ap(1)
       ape(38)=ap(13)
       ape(39)=ap(13)

       DO k = 1,13
          i = (k * 3) - 1
          ape(i) = ap(k)
          END DO

       DO k = 1,12
          i = k * 3
          ape(i) = (ap(k)*2 + ap(k+1))/3.0
          END DO

       DO k = 2,13
          i = (k * 3) - 2
          ape(i) = (ap(k-1) + ap(k)*2)/3.0
          END DO

C     FROM 3-HOURLY TO HOURLY ap (old version without interpolation)
c      i = 1
c      DO 10 k = 1,13
c         DO j = 1,3
c            ape(i) = ap(k)
c            i = i + 1
c            END DO
c10    CONTINUE

C     TO OBTAIN THE INTEGRAL OF ap.
C     INTEGRAL OF ap

      if(ut.eq.24) ut=0
      IF (ut .EQ. 0 .OR. ut .EQ. 3 .OR. ut .EQ. 6 .OR. ut .EQ. 9 .OR.
     1ut .EQ. 12 .OR. ut .EQ. 15 .OR. ut .EQ. 18 .OR. ut .EQ. 21) THEN
          k = 1
      ELSE IF (ut .EQ. 1 .OR. ut .EQ. 4 .OR. ut .EQ. 7 .OR. ut .EQ. 10
     1.OR.ut .EQ. 13 .OR. ut .EQ. 16 .OR. ut .EQ. 19 .OR. ut .EQ. 22)
     2THEN
          k = 2
      ELSE IF (ut .EQ. 2 .OR. ut .EQ. 5 .OR. ut .EQ. 8 .OR. ut .EQ. 11
     1.OR. ut .EQ. 14 .OR. ut .EQ. 17 .OR. ut .EQ. 20 .OR. ut .EQ. 23)
     2THEN
          k = 3

      ELSE

          WRITE (6,*)' '
          WRITE (6,*)' '
          WRITE (6,*)'  Wrong Universal Time value -------- >>', ut
          WRITE (6,*)' '
          GOTO 100

      END IF

      rap = 0

      DO j = 1,36
      rap = rap + fap(j) * ape(k+j)
      END DO

      if(rap.le.200.)then
      cf=1.0
      goto 100
      end if

      if(doy.gt.366.or.doy.lt.1)then
          WRITE (6,*)' '
          WRITE (6,*)' '
          WRITE (6,*)' '
          WRITE (6,*)'      Wrong Day of Year value --- >>', doy
          WRITE (6,*)' '
          GOTO 100
      end if

      if(rgma.gt.90.0.or.rgma.lt.-90.0)then
          WRITE (6,*)' '
          WRITE (6,*)' '
          WRITE (6,*)' '
          WRITE (6,*)'   Wrong GEOMAGNETIC LATITUDE value --- >>', rgma
          WRITE (6,*)' '
          GOTO 100
      end if

c      write(6,*)rgma

      dayno=doy
      if(rgma.lt.0.0)then
      dayno=doy+172
      if(dayno.gt.365)dayno=dayno-365
      end if

      if (dayno.ge.82) rs=(dayno-82.)/45.6+1.
      if (dayno.lt.82) rs=(dayno+283.)/45.6+1.
      s1=rs
      facs=rs-s1
      s2=s1+1
      if(s2.eq.9) s2=1
c      write(6,*)s1,s2,rs

      rgma = abs(rgma)

      rl=(rgma+10.)/20.+1
      if(rl.eq.6.0)rl=5.9
      l1=rl
      facl=rl-l1
      l2=l1+1
c      write(6,*)l1,l2,rl

C     FACTORS CALCULATIONS

      if(rap.lt.300.)then
      rapf=300.
      n1=code(s1,l1)
      cf1=c4(n1)*(rapf**4)+c3(n1) * (rapf**3) + c2(n1) * (rapf**2) +
     1c1(n1) * rapf + c0(n1)
      n2=code(s1,l2)
      cf2=c4(n2)*(rapf**4)+c3(n2) * (rapf**3) + c2(n2) * (rapf**2) +
     1c1(n2) * rapf + c0(n2)
      n3=code(s2,l1)
      cf3=c4(n3)*(rapf**4)+c3(n3) * (rapf**3) + c2(n3) * (rapf**2) +
     1c1(n3) * rapf + c0(n3)
      n4=code(s2,l2)
      cf4=c4(n4)*(rapf**4)+c3(n4) * (rapf**3) + c2(n4) * (rapf**2) +
     1c1(n4) * rapf + c0(n4)

C     INTERPOLATION

      cf300=cf1*(1 - facs) * (1 - facl) + cf2 * (1 - facs) * (facl) +
     *cf3 * (facs) * (1 - facl) + cf4 * (facs) * (facl)

      cf = (cf300-1.0)*rap/100.-2.*cf300+3.
      goto 100
      end if

      n1=code(s1,l1)
c      write(6,*)n1
      cf1 = c4(n1) * (rap**4) + c3(n1) * (rap**3) + c2(n1) * (rap**2) +
     1c1(n1) * rap + c0(n1)
      n2=code(s1,l2)
      cf2 = c4(n2) * (rap**4) + c3(n2) * (rap**3) + c2(n2) * (rap**2) +
     1c1(n2) * rap + c0(n2)
      n3=code(s2,l1)
      cf3 = c4(n3) * (rap**4) + c3(n3) * (rap**3) + c2(n3) * (rap**2) +
     1c1(n3) * rap + c0(n3)
      n4=code(s2,l2)
      cf4 = c4(n4) * (rap**4) + c3(n4) * (rap**3) + c2(n4) * (rap**2) +
     1c1(n4) * rap + c0(n4)

C     INTERPOLATION

      cf = cf1 * (1 - facs) * (1 - facl) + cf2 * (1 - facs) * (facl) +
     *cf3 * (facs) * (1 - facl) + cf4 * (facs) * (facl)

100   CONTINUE

      RETURN

      END
C
c
c
      subroutine igrf_dip(xlat,xlong,year,height,dec,dip,dipl,ymodip)
c-----------------------------------------------------------------------        
c INPUT:
c    xlat      geodatic latitude in degrees
c    xlong     geodatic longitude in degrees
c    year      decimal year (year+month/12.0-0.5 or 
c                  year+day-of-year/365 or ../366 if leap year) 
c    height    height in km
c OUTPUT:
c    dec       magnetic declination in degrees
c    dip       magnetic inclination (dip) in degrees
c    dipl      dip latitude in degrees
c    ymodip    modified dip latitude = asin{dip/sqrt[dip^2+cos(LATI)]} 
c-----------------------------------------------------------------------        

      COMMON/IGRF1/     UMR,ERA,AQUAD,BQUAD
C
      CALL INITIZE
c
C----------------CALCULATE PROFILES-----------------------------------
c
		xlati = xlat
		xlongi = xlong
		h = height
        CALL FELDCOF(YEAR,DIMO)
        CALL FELDG(XLATI,XLONGI,H,BNORTH,BEAST,BDOWN,BABS)
          DECARG=BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH)
          IF(ABS(DECARG).GT.1.) DECARG=SIGN(1.,DECARG)
        DEC=ASIN(DECARG)
          BDBA=BDOWN/BABS
          IF(ABS(BDBA).GT.1.) BDBA=SIGN(1.,BDBA)
        DIP=ASIN(BDBA)
          dipdiv=DIP/SQRT(DIP*DIP+cos(XLATI*UMR))
          IF(ABS(dipdiv).GT.1.) dipdiv=SIGN(1.,dipdiv)
        SMODIP=ASIN(dipdiv)
        
c       DIPL1=ATAN(0.5*TAN(DIP))/UMR

      DIPL=ATAN(BDOWN/2.0/sqrt(BNORTH*BNORTH+BEAST*BEAST))/umr
      YMODIP=SMODIP/UMR                            
      DEC=DEC/UMR
      DIP=DIP/UMR
      RETURN
      END
c
c
        SUBROUTINE FELDCOF(YEAR,DIMO)
c-----------------------------------------------------------------------        
C  DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
C
C       INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
C                       BE CALCULATED
C       OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED 
C                       TO EARTH'S RADIUS) AT THE TIME (YEAR)
C 05/31/2000 updated to IGRF-2000 version (###) 
C 03/24/2000 updated to IGRF-2005 version (###) 
C 07/22/2009 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
C 02/26/2010 update to IGRF-11 (2010) (###)  
C 10/05/2011 added COMMON/DIPOL/ for MLT computation in DPMTRX (IRIFUN)
C 02/10/2015 update to IGRF-12 (2015) (###)
c-----------------------------------------------------------------------        
        CHARACTER*13    FILMOD, FIL1, FIL2           
C ### FILMOD, DTEMOD array-size is number of IGRF maps
        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(16)
        DIMENSION		DTEMOD(16)
        DOUBLE PRECISION X,F0,F 
        COMMON/MODEL/   NMAX,TIME,GH1,FIL1
        COMMON/IGRF1/   UMR,ERAD,AQUAD,BQUAD
        COMMON/DIPOL/	GHI1,GHI2,GHI3
C ### updated coefficient file names and corresponding years
        DATA  FILMOD   / 'dgrf1945.dat','dgrf1950.dat','dgrf1955.dat',           
     1    'dgrf1960.dat','dgrf1965.dat','dgrf1970.dat','dgrf1975.dat',
     2    'dgrf1980.dat','dgrf1985.dat','dgrf1990.dat','dgrf1995.dat',
     3    'dgrf2000.dat','dgrf2005.dat','dgrf2010.dat','igrf2015.dat',
     4    'igrf2015s.dat'/
        DATA  DTEMOD / 1945., 1950., 1955., 1960., 1965.,           
     1   1970., 1975., 1980., 1985., 1990., 1995., 2000.,2005.,
     2   2010., 2015., 2020./      
C
C ### numye is number of IGRF coefficient files minus 1
C
        NUMYE=15
C
C  IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
C  IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
C
        IU = 14
        IS = 0
C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR
        TIME = YEAR
        IYEA = INT(YEAR/5.)*5
        L = (IYEA - 1945)/5 + 1
        IF(L.LT.1) L=1
        IF(L.GT.NUMYE) L=NUMYE         
        DTE1 = DTEMOD(L)   
        FIL1 = FILMOD(L)   
        DTE2 = DTEMOD(L+1) 
        FIL2 = FILMOD(L+1) 
C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
            IF (IER .NE. 0) STOP                           
        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
            IF (IER .NE. 0) STOP
C-- DETERMINE IGRF COEFFICIENTS FOR YEAR
        IF (L .LE. NUMYE-1) THEN                        
          CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, 
     1          NMAX2, GH2, NMAX, GHA)                        
        ELSE               
          CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,     
     1          GH2, NMAX, GHA)                                    
        ENDIF 
C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
        F0=0.D0
        DO 1234 J=1,3
           F = GHA(J) * 1.D-5
           F0 = F0 + F * F
1234    CONTINUE
        DIMO = DSQRT(F0)
        GHI1=GHA(1)                                         
        GHI2=GHA(2)                                         
        GHI3=GHA(3)

        GH1(1) =  0.0
        I=2          
        F0=1.D-5                
        IF(IS.EQ.0) F0=-F0 
        SQRT2=SQRT(2.)      

      DO 9 N=1,NMAX           
        X = N
        F0 = F0 * X * X / (4.D0 * X - 2.D0)               
        IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
        F = F0 * 0.5D0                                    
        IF(IS.EQ.0) F = F * SQRT2
        GH1(I) = GHA(I-1) * F0
        I = I+1                                         
      DO 9 M=1,N                                    
        F = F * (X + M) / (X - M + 1.D0)                 
        IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))             
        GH1(I) = GHA(I-1) * F
        GH1(I+1) = GHA(I) * F
        I=I+2
9     CONTINUE 
        RETURN
        END
C
C
        SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                                                                                           
C ===============================================================               
C       Reads spherical harmonic coefficients from the specified     
C       file into an array.                                          
C       Input:                                                       
C           IU    - Logical unit number                              
C           FSPEC - File specification                               
C       Output:                                                      
C           NMAX  - Maximum degree and order of model                
C           ERAD  - Earth's radius associated with the spherical     
C                   harmonic coefficients, in the same units as      
C                   elevation                                        
C           GH    - Schmidt quasi-normal internal spherical          
C                   harmonic coefficients                            
C           IER   - Error number: =  0, no error                     
C                                 = -2, records out of order         
C                                 = FORTRAN run-time error number    
C ===============================================================               
                                                                                
        CHARACTER  FSPEC*(*), FOUT*80                                    
        DIMENSION       GH(196)
        LOGICAL		mess 
        COMMON/iounit/konsol,mess        
        do 1 j=1,196  
1          GH(j)=0.0

C ---------------------------------------------------------------               
C       Open coefficient file. Read past first header record.        
C       Read degree and order of model and Earth's radius.           
C ---------------------------------------------------------------               
        WRITE(FOUT,667) FSPEC
 667    FORMAT(A13)
c 667    FORMAT('/var/www/omniweb/cgi/vitmo/IRI/',A13)
        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)     
        READ (IU, *, IOSTAT=IER, ERR=999)                            
        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD, XMYEAR 
        nm=nmax*(nmax+2)                
        READ (IU, *, IOSTAT=IER, ERR=999) (GH(i),i=1,nm) 
        goto 888 
               
999     write(konsol,100) FOUT
100     FORMAT('Error while reading ',A13)

888     CLOSE (IU)                                                                                                                                   
        RETURN                                                       
        END                                                          
C
C
        SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,          
     1                        NMAX2, GH2, NMAX, GH)                  
                                                                                
C ===============================================================               
C                                                                               
C       Version 1.01                                                 
C                                                                               
C       Interpolates linearly, in time, between two spherical        
C       harmonic models.                                             
C                                                                               
C       Input:                                                       
C           DATE  - Date of resulting model (in decimal year)        
C           DTE1  - Date of earlier model                            
C           NMAX1 - Maximum degree and order of earlier model        
C           GH1   - Schmidt quasi-normal internal spherical          
C                   harmonic coefficients of earlier model           
C           DTE2  - Date of later model                              
C           NMAX2 - Maximum degree and order of later model          
C           GH2   - Schmidt quasi-normal internal spherical          
C                   harmonic coefficients of later model             
C                                                                               
C       Output:                                                      
C           GH    - Coefficients of resulting model                  
C           NMAX  - Maximum degree and order of resulting model      
C                                                                               
C       A. Zunde                                                     
C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
C                                                                               
C ===============================================================               
                                                                                
        DIMENSION       GH1(*), GH2(*), GH(*)                        
                                                                                
C ---------------------------------------------------------------               
C       The coefficients (GH) of the resulting model, at date        
C       DATE, are computed by linearly interpolating between the     
C       coefficients of the earlier model (GH1), at date DTE1,       
C       and those of the later model (GH2), at date DTE2. If one     
C       model is smaller than the other, the interpolation is        
C       performed with the missing coefficients assumed to be 0.     
C ---------------------------------------------------------------               
                                                                                
        FACTOR = (DATE - DTE1) / (DTE2 - DTE1)                       
                                                                                
        IF (NMAX1 .EQ. NMAX2) THEN                                   
            K = NMAX1 * (NMAX1 + 2)                                  
            NMAX = NMAX1                                             
        ELSE IF (NMAX1 .GT. NMAX2) THEN                              
            K = NMAX2 * (NMAX2 + 2)                                  
            L = NMAX1 * (NMAX1 + 2)                                  
            DO 1122 I = K + 1, L                                          
1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
            NMAX = NMAX1                                             
        ELSE                                                         
            K = NMAX1 * (NMAX1 + 2)                                  
            L = NMAX2 * (NMAX2 + 2)                                  
            DO 1133 I = K + 1, L                                          
1133            GH(I) = FACTOR * GH2(I)                              
            NMAX = NMAX2                                             
        ENDIF                                                        
                                                                                
        DO 1144 I = 1, K                                                  
1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
                                                                                
        RETURN                                                       
        END                                                          
C
C
        SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,           
     1                        GH2, NMAX, GH)                           
                                                                                
C ===============================================================               
C                                                                               
C       Version 1.01                                                   
C                                                                               
C       Extrapolates linearly a spherical harmonic model with a        
C       rate-of-change model.                                          
C                                                                               
C       Input:                                                         
C           DATE  - Date of resulting model (in decimal year)          
C           DTE1  - Date of base model                                 
C           NMAX1 - Maximum degree and order of base model             
C           GH1   - Schmidt quasi-normal internal spherical            
C                   harmonic coefficients of base model                
C           NMAX2 - Maximum degree and order of rate-of-change         
C                   model                                              
C           GH2   - Schmidt quasi-normal internal spherical            
C                   harmonic coefficients of rate-of-change model      
C                                                                               
C       Output:                                                        
C           GH    - Coefficients of resulting model                    
C           NMAX  - Maximum degree and order of resulting model        
C                                                                               
C       A. Zunde                                                       
C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      
C                                                                               
C ===============================================================               
                                                                                
        DIMENSION       GH1(*), GH2(*), GH(*)                        
                                                                                
C ---------------------------------------------------------------               
C       The coefficients (GH) of the resulting model, at date          
C       DATE, are computed by linearly extrapolating the coef-         
C       ficients of the base model (GH1), at date DTE1, using          
C       those of the rate-of-change model (GH2), at date DTE2. If      
C       one model is smaller than the other, the extrapolation is      
C       performed with the missing coefficients assumed to be 0.       
C ---------------------------------------------------------------               
                                                                                
        FACTOR = (DATE - DTE1)                                         
                                                                                
        IF (NMAX1 .EQ. NMAX2) THEN                                     
            K = NMAX1 * (NMAX1 + 2)                                    
            NMAX = NMAX1                                               
        ELSE IF (NMAX1 .GT. NMAX2) THEN                                
            K = NMAX2 * (NMAX2 + 2)                                    
            L = NMAX1 * (NMAX1 + 2)                                    
            DO 1155 I = K + 1, L                                            
1155            GH(I) = GH1(I)                                         
            NMAX = NMAX1                                               
        ELSE                                                           
            K = NMAX1 * (NMAX1 + 2)                                    
            L = NMAX2 * (NMAX2 + 2)                                    
            DO 1166 I = K + 1, L                                            
1166            GH(I) = FACTOR * GH2(I)                                
            NMAX = NMAX2                                               
        ENDIF                                                          
                                                                                
        DO 1177 I = 1, K                                                    
1177        GH(I) = GH1(I) + FACTOR * GH2(I)                           
                                                                                
        RETURN                                                         
        END                                                            
C
C
        SUBROUTINE INITIZE
C----------------------------------------------------------------
C Initializes the parameters in COMMON/IGRF1/
C
C       UMR     = ATAN(1.0)*4./180.   <DEGREE>*UMR=<RADIANT>
C       ERA     EARTH RADIUS FOR NORMALIZATION OF CARTESIAN 
C                       COORDINATES (6371.2 KM) 
C       EREQU   MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
C       ERPOL   MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
C       AQUAD   SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
C       BQUAD   SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
C
C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL 
C ASTRONOMICAL UNION .
C-----------------------------------------------------------------
        COMMON/IGRF1/   UMR,ERA,AQUAD,BQUAD
        ERA=6371.2
        EREQU=6378.16
        ERPOL=6356.775
        AQUAD=EREQU*EREQU
        BQUAD=ERPOL*ERPOL
        UMR=ATAN(1.0)*4./180.
        RETURN
        END
C
C
      SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)           
c-----------------------------------------------------------------------        
C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61, 
C      1970.
c-----------------------------------------------------------------------        
C CHANGES (D. BILITZA, NOV 87):
C   - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA
C   - CALCULATES DIPOL MOMENT
C 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
c-----------------------------------------------------------------------        
C  INPUT:  ENTRY POINT FELDG
C               GLAT  GEODETIC LATITUDE IN DEGREES (NORTH)
C               GLON  GEODETIC LONGITUDE IN DEGREES (EAST)
C               ALT   ALTITUDE IN KM ABOVE SEA LEVEL
C
C          ENTRY POINT FELDC
C               V(3)  CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
C                       X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
C                       Y-AXIS POINTING TO EQUATOR AT 90 LONG.
C                       Z-AXIS POINTING TO NORTH POLE
C
C          COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED
C            IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG.
C       
C          COMMON /MODEL/ AND /IGRF1/
C               UMR     = ATAN(1.0)*4./180.   <DEGREE>*UMR=<RADIANT>
C               ERA     EARTH RADIUS FOR NORMALIZATION OF CARTESIAN 
C                       COORDINATES (6371.2 KM)
C               AQUAD, BQUAD   SQUARE OF MAJOR AND MINOR HALF AXIS OF 
C                       EARTH ELLIPSOID AS RECOMMENDED BY INTERNAT. 
C                       ASTRONOMICAL UNION (6378.160, 6356.775 KM).
C               NMAX    MAXIMUM ORDER OF SPHERICAL HARMONICS
C               TIME    YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC 
C                       FIELD IS TO BE CALCULATED
C               G(M)    NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF)
C                       M=NMAX*(NMAX+2)
c-----------------------------------------------------------------------        
C  OUTPUT: BABS   MAGNETIC FIELD STRENGTH IN GAUSS
C          BNORTH, BEAST, BDOWN   COMPONENTS OF THE FIELD WITH RESPECT
C                 TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS
C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
C                 AND DOWNWARD.   
C-----------------------------------------------------------------------
      DIMENSION         V(3),B(3)   
      CHARACTER*13      NAME
      COMMON/IGRF2/	XI(3),H(196)
      COMMON/MODEL/     NMAX,TIME,G(196),NAME  
      COMMON/IGRF1/     UMR,ERA,AQUAD,BQUAD

C
C-- IS RECORDS ENTRY POINT
C
C*****ENTRY POINT  FELDG  TO BE USED WITH GEODETIC CO-ORDINATES         
      IS=1                                                              
      RLAT=GLAT*UMR
      CT=SIN(RLAT)                                                      
      ST=COS(RLAT)                                                      
      D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)                                 
      RLON=GLON*UMR
      CP=COS(RLON)                                                      
      SP=SIN(RLON)                                                      
      ZZZ=(ALT+BQUAD/D)*CT/ERA
      RHO=(ALT+AQUAD/D)*ST/ERA
      XXX=RHO*CP                                                       
      YYY=RHO*SP                                                       
      GOTO 10                                                            

C*****ENTRY POINT  FELDC  TO BE USED WITH CARTESIAN CO-ORDINATES        
      ENTRY FELDC(V,B)                                                  
      IS=2                                                              
      XXX=V(1)                                                          
      YYY=V(2)                                                          
      ZZZ=V(3)                                                          
10    RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ) 
      XI(1)=XXX*RQ                                                      
      XI(2)=YYY*RQ                                                      
      XI(3)=ZZZ*RQ                                                      
      GOTO 20                                                            

C*****ENTRY POINT  FELDI  USED FOR L COMPUTATION                        
      ENTRY FELDI                                                       
      IS=3                                                              
20    IHMAX=NMAX*NMAX+1                                                 
      LAST=IHMAX+NMAX+NMAX                                              
      IMAX=NMAX+NMAX-1                                                  
      DO 8 I=IHMAX,LAST                                                 
8     	H(I)=G(I)                                                         
      DO 6 K=1,3,2                                                      
      	I=IMAX                                                            
      	IH=IHMAX                                                          
1     	IL=IH-I                                                           
      	F=2./FLOAT(I-K+2)                                                 
      	X=XI(1)*F                                                         
      	Y=XI(2)*F                                                         
      	Z=XI(3)*(F+F)                                                     
      	I=I-2                                                             
      	IF(I-1) 5,4,2
      	                                                      
2     	DO 3 M=3,I,2                                                      
      		H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-          
     A                H(IH+M-1))-Y*(H(IH+M+2)+H(IH+M-2))           
3     		H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)-                
     A                H(IH+M-2))+Y*(H(IH+M+3)+H(IH+M-1))
                      
4     	H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH))             
      	H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH))             
5     	H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2))                      
      	IH=IL                                                             
      	IF(I.GE.K) GOTO 1                                                   
6     	CONTINUE
                                                          
      IF(IS.EQ.3) RETURN 
                                                      
      S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))                   
      T=(RQ+RQ)*SQRT(RQ)                                                
      BXXX=T*(H(3)-S*XXX)                                               
      BYYY=T*(H(4)-S*YYY)                                               
      BZZZ=T*(H(2)-S*ZZZ) 
                                                    
      IF(IS.EQ.2) GOTO 7                                                  
      BABS=SQRT(BXXX*BXXX+BYYY*BYYY+BZZZ*BZZZ)
      BEAST=BYYY*CP-BXXX*SP                                             
      BRHO=BYYY*SP+BXXX*CP                                              
      BNORTH=BZZZ*ST-BRHO*CT                                            
      BDOWN=-BZZZ*CT-BRHO*ST                                            
      RETURN 
                                                                 
7     B(1)=BXXX                                                         
      B(2)=BYYY                                                         
      B(3)=BZZZ                                                         
      RETURN                                                            
      END                                                               
C
C
